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1° Introduction and Goals 


These notes are intended to give computer graphics programmers 
and artists an introduction to methods of simulating, animating, and 
rendering ocean water environments. CG water has become a com- 
‘mon tool in visual effects work at all levels of computer graphics, 
{rom print media to feature films. Several commercial products are 
available for nearly any computer platform and work environment. 
A few visual effects companies continue to extend and improve 
these tools, seeking to generate higher quality surface geometry, 
complex inteatctions, and more compelling imagery. In order for 
an artist to exploit these tools to maximum benefit, itis important 
that he or she become familiar with concepts, terminology, a litle 
oceanography, and the present state ofthe at 

‘As demonstrated by the pioneering efforts in the films Water- 
World and Titanic, as well as several other films made since about 
1995, images of eg water can be generated with a high degree of 
realism. However, this level of realism has been mostly limited 
to relatively calm, nice ocean conditions. Conditions with large 
amounts of spray, breaking waves, foam, splashing, and wakes are 
improving and approaching the same realistic look 

Three general approaches are currently popular in computer 
graphics for simulating fluid motion, and water surfaces in partic- 
lular. The three methods are all related in some way to the basic 
Navier-Stokes equations at the heart of many applications, Com- 
putational Fluid Dynamics (CFD) is one of the methods which has 
recieved a great amount of attention lately, While many versions of| 
this technology have been in existence for some time, recent papers 
by Stam (0), Fedkiw and Foster (2), and many others have demon- 
strated thatthe numerical computation discipline and the computer 
graphics discipline have connected well enough to produce useful, 
interesting, and sometimes beautiful results. The primary draw- 
back of CFD methods is thatthe computations are performed on 
data structured in a 2D or 3D grid, called an Eulerian framework. 
This gridded architecture limits the combined extent and flow detail 
that can be computed. At present for example, itis practical to sim- 
ulate with CFD waves breaking as they approach a shallow beach, 
However, that simulation is not able to simulate a bay-sized region 
of ocean while also simulating the breaking down to the detail of| 
spray formation, simply because the number of grid points needed 
is prohibitively large, 

‘A second method that is shows great promise is called Smoothed 
Particle Hydrodynamics (SPH). This is a completely different ap- 
proach to solving the Navier-Stokes equations. SPH imagines that 
the volume ofa fluid is composed of small overlapping regions, the 
center of each region carrying some amount of mass and momen- 
tum, The regions are allowed to move about within the uid and ex- 
perience forces due to pressure, strain, gravity, and others. But now 
the center of each region acts as a particle, and the Navier-Stokes 
equations are converted into equations of motion for discrete par- 
ticles. This approach is called a Lagragian framework (as opposed 
to the Eulerian framework in CFD). The fluid volume is not bound 
to grid geometry. SPH is a very useful method of simulation for 
situations in which there is significant splashing or explosions, and 
ta evtn bets asad to slated racking i ste E. Uaing fn 
plicit methods to construct surface forthe uid, standard computer 
graphics applications such as pouring water have been achieved (3) 

Finally, the third method is the one that is the focus of these 
notes, and unlike the CFD and SPH methods, this one is focused 
6 the more narrow goal of simulating the motion of the surface 
of a body of water, Surface simulations are commonly generated 
from the CFD and SPH methods, but the surface is generated by 
algorithms that are added to those fluid simulations, for example by 
tracking a type of implicit surface called a level set. In choosing 
to focus on the surface structure and motion, we eliminate much of 
the computation and resolution limits in the CFD and SPH meth- 


‘ods. Of course when we eliminate those computations we have lost 
certain types of realistic motion, most notably the breakup of the 
surface with strong changes in topology. In place of that compu- 
lation we substitute a mixture of knowledge of oceanographic phe- 
nomenology and computational flexibility to achieve realistic types 
‘of surfaces that cannot practically be achieved by the others. In that 
sense, this phemonological approach should be considered comple~ 
mentary ~ not competitive ~ to the CFD and SPH methods, since 
all three work best in different regimes of fluid motion. 
Broadly, the reader should come away from this material with, 


1. an understanding of the important physical concepts for ocean 
surface propagation, most notably the concept of dispersion 
and types of dispersion relationships. 


aan understanding of some algorithms that generate/animate 
water surface height fields suitable for modeling waves as big 
as storm surges and as small as tiny capillaries; 


3. an understanding of the basic optical processes of reflection 
and refraction from a water surface; 


4. an introduction to the color filtering behavior of ocean water; 


5. an introduction to complex lighting effects known as caus- 
tics and godrays, produced when sunlight passes through the 
rough surface into the water Volume underneath; and 


6. some rules of thumb for which choices make nice looking 
images and what are the tradeoffs of quality versus compu- 
tational resources. Some example shaders are provided, and 
‘example renderings demonstrate the content of the discussion, 


Before diving into it, I frst want to be more concrete about what, 
aspect of the ocean environment we cover (or not cover) in these 
notes. Figure[T]}is a rendering of an oceanseape produced from mod- 
cls of water, ai, and clouds. Light from the cloues is reflected from 
the surface. On the extreme left, sun gliter is also present. The 
‘generally bluish color of the water is due to the reflection of blue 
skylight, and to light coming out of the water after scattering from 
the volume. Although these notes do not tackle the modeling and 
rendering of clouds and air, there is a discussion of how skylight 
from the clouds and air is reflected from, or refracted through, the 
water surface. These notes will tell you how to make a height-field 
displacement-mapped surface for the ocean waves with the detail 
and quality shown in the figure. The notes also discuss several ef- 
fects of the underwater environment and how to model/render them. 
‘The primary four effects are sunbeams (also called godrays), caus- 
tics on underwater surfaces, blurring by the scattering of light, and 
color filtering. 

‘There are also many other complex and interesting aspects ofthe 
‘ocean environment that will not be covered. These include break- 
ing waves, spray, foam, wakes around objects in the water. splashes 
from bodies that impact the surface. and global illumination of the 
entire ocean-atmosphere environment. There is substantial research 
‘underway on these topics, and so it is possible that future versions 
of this or other lecture notes will include them. I have included a 
brief section on advanced modifications to the basic wave height al- 
{gorithm that produce choppy waves. The modification could feasi 
bly lead toa complete description ofthe surface portion of breaking 
waves, and possibly serve to drive the spray and foam dynamics as 
well 

‘There is, of course, a substantial body of literature on ocean sur- 
face simulation and animation, both in computer graphics circles 
and in oceanography. One ofthe first descriptions of water wavesin 
computer graphics was by Fournier and Reeves(T2} . who modeled 
a shoreline with waves coming up on it using a water surface model 
called Gerstner waves. In that same issue, Darwin Peachey(13) 


2. RADIOSITY OF THE OCEAN ENVIRONMENT 


Figure 1: Rendered image of an oceanscape. 


presented a variation on this approach using basis shapes other than 
sinusoids. 

In the oceanographic literature, ocean optics became an inten- 
sive topic of research in the 1940s. $.Q. Duntley published{T7} in 
1963 papers containing optical data of relevance to computer graph- 
ies, Work continues today. The field of optical oceanography has 
grown into a mature quantitative science with subdisciplines and 
many different applications. One excellent review of the state of 
the science was written by Curtis Mobley() 

In these lectures the approach we take to creating surface waves, 
is close to the one outlined by Masten, Watterberg, and Mareda(T), 
although the technique had been in use for many years prior to 
their paper in the optical oceanography community. This approach 
synthesizes a patch of ocean waves from a Fast Fourier Transform 
(FFT) prescription, with user-controllable size and resolution, and 
which can be tiled seamlessly over a larger domain. The patch con- 
tains many octaves of sinusoidal waves that all add up at each point 
to produce the synthesized height. The mixture of sinusoidal am- 
plitudes and phases however, comes from statistical, emperically- 
based models of the ocean. What makes these sinusoids look like 
‘waves and not just a bunch of sine waves is the large collection of 
sinusoids that are used, the relative amplitudes of the sinusoids, and 
their animation using the dispersion relation. We examine the im- 
pact of the number of sinusoids and resolution on the quality of the 
rendered image. 

In the next section we begin the discussion ofthe ocean environ- 
ment with a broad introduction to the global illamination problem, 
The radiosity equations for this environment look much like those 
of any other radiosity problem, although the volumetric character of 
some of the environmental components complicate a general imple- 
mentation considerably. However, we simplify the issues by ignor- 
ing some interactions and replacing others with models generated 
by remote sensing data 

Practical methods are presented in section[I}for creating realiza- 
tions of ocean surfaces. We present two metiods, one based on a 
simple model of water structure and movement, and one based on 
summing up large numbers of sine waves with amplitudes that are 
related to each other based on experimental evidence. This sec- 
cond method carries out the sum using the technique of Fast Fourier 
‘Transformation (ft), and has been used effectively in projects for 
commercials, television, and motion pictures. 

‘Afler the discussion of the structure and animation of the water 
surface, we focus on the optical properties of water relevant to the 


‘graphics problem. First, we discuss the interaction at the air-water 
interface: reflection and refraction, This leaves us with a simple 
but effective Renderman-style shader suitable for rendering water 
surfaces in BMRT, for example. Next, the optical characteristics of 
the underwater environment are explored. 

Finally, please remember that these notes are a living 
document, Some of the discussion of the various top- 
ies is still very limited and incomplete. If you find a 
problem or have additional questions, please feel free to 


contact me at The 


latest version of course documents are hosted at 


ip: 7heww.finclighivisuallechnology.com| 


2  Radiosity of the Ocean Environment 


‘The ocean environment, for our purposes, consists of only four 
components: The water surface, the air, the sun, and the water 
yolume below the surface. In this section we trace the flow of 
light through the environment, both mathematically and schemat- 
ically, from the light source to the camera. In general, the radiosity 
‘equations here are as coupled as any other radiosity problem. To a 
reasonable degree, however, the coupling can be truncated and the 
simplified radiosity problem has a relatively fast solution, 

‘The light seen by a camera is dependent on the flow of light en- 
ergy from the source(s) (ie. the sun and sky) to the surface and 
into the camera. In addition to specular reflection of direct sun- 
light and skylight from the surface, some fraction of the incident 
light is transmitted through the surface. Ultimately, a fraction of the 
transmitted light is scattered by the water volume back up through 
the interface and into the air. Some of the light that is reflected or 
refracted at the surface may strike the surface a second time, pro- 
ducing more reflection and refraction events. Under some viewing 
conditions, multiple reflections and refractions can have a notice- 
able impact on images. For our part however, we will ignore more 
than one reflection or refraction from the surface at a time. This 
not only makes the algorithms and computation easier and faster, 
but also is reasonably accurate in most viewing conditions and pro- 
duces visually realistic imagery, 

At any point in the environment above the surface, including at 
the camera, the total light intensity (radiance) coming from any di 
rection has three contributions: 


Lanove =rhs+rLa + tule , a 
with the following definitions of the terms 


ris the Fresnel reflectivity for reflection from a spot on the surface 
of the ocean to the camera, 


fy is the transmission coefficient forthe light Lry coming up from 
the ocean volume, refracted at the surface into the camera, 


Ls is the amount of light coming directly from the sun, through 
the atmosphere, to the spot on the ocean surface where it is 
reflected by the surface to the camera. 


Lis the (diffuse) atmospheric skylight 


Ly is the light just below the surface that is transmitted through 
the surface into the air. 


Equation[T]has intentionally been written in a shorthand way that 
hhides the dependences on position in space and the direction the 
light is traveling. 

While equation [T]appears to have a relatively simple structure, 
the terms Ls, a, and Lv ean in principle have complex depen- 
dencies on each other, as well on the reflectivity and transmissivity. 
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There is a large body of research literature investigating these de- 
pendencies in detail (19), but we will not at this point pursue these 
quantitative methods. But we can elaborate further on the coupling 
while continuing with the same shorthand notation, The direct light 
from the sun Ls is 


Ls = Lroaexp{-r} 2 


where Lo. is the intensity of the direct sunlight at the top of the 
atmosphere, and r is the “optical thickness” of the atmosphere for 
the direction of the sunlight and the point on the earth. Both the 
diffuse atmospheric skylight L. and the upwelling light Ly can be 
‘written as the sum of two terms: 


La = La(Ls)+£4(Lu) cc) 
Le = Lb(Ls)+Li(La) 4) 


These equations reveal the potential complexity of the problem. 
While both £4 and Ly depend on the direct sunlight. they also 
depend on each other. For example, the total amount of light pene- 
trating into the ocean comes from the direct sunlight and from the 
atmospheric sunlight. Some of the light coming into the ocean is 
scattered by particulates and molecules in the ocean, back up into 
the atmosphere, Some of that upwelling light in turn is scattered in 
the atmosphere and becomes a part of the skylight shining on the 
surface, and on and on. This is a classic problem in radiosity. It 
is not particularly special for this case, as opposed to other radios- 
ity problems, except perhaps for the fact that the upwelling light 
is difficult to compute because it comes from volumetric multiple 
scattering, 

‘Our approach, for the purposes of these notes, to solving this 
radiosity problem is straightforward: take the skylight to depend 
only on the light from the sun, since the upwelling contribution 
represents a “tertiary” dependence on the sunlight; and completely 
replace the equation for Ly with an empirical formula, based on 
scientific observations of the oceans, that depends only on the di- 
rect sunlight and a few other parameters that dictate water type and 
clarity. 

Under the water surface, the radiosity equation has the schematic 
form 

Laevow = thn +tLr + Les + Lar 6) 


with the meaning 


{is the Fresnel transmissivity for transmission through the water 
surface at each point and angle on the surface. 


Lp The “direct” light from the sun that penetrates into the water. 


Li The “indirect” light from the atmosphere that penetrates into 
the water. 


Ls The single-scattered light, from both the sun and the atmo- 
sphere, that is scattered once in the water volume before ar- 
riving at any point 


Lay The multiply-scattered light. This is the single-scattered light 
that undergoes more scattering events in the volume. 


Just as for the above water case, these terms are all related to each 
other is relative complex ways. For example, the single scattered 
light depends on the direct and indirect light: 


Lss = P(tLs) + P(tLo) © 
with the quantity P being a linear functional operator of its argu- 


‘ment, containing information about the single scattering event and 
the aitenuation of the scattered light as it passes from the scatter 
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Figure 2: Illustration of multiple reflections and transmission 
through the air-water interface. 


point to the camera. Similarly, the multiply-scattered light is de- 
pendent on the single scattered: 


Lar = G(tLs) + G(tLo) a) 


‘The functional schematic quantities P and G are related, since mul- 
tiple scattering is just a series of single scatters. Formally, the two 
hhave an operator dependence that has the form 


~ P@P®@exp(P) (8) 


[At this point, the schematic representation may have outlived its 
tusefullness because of the complex (and not here defined) meaning 
of the convolution-like operator ©, and because the expression for 
G in terms of P has created an even more schematic view in terms 
of an exponentisted P. So for now we will leave the schematic 
representation, and journey on with more concrete quantities the 
rest of the way through 

‘The formal schematic discussion put forward here does have a 
mathematically and physically precise counterpart. The field of 
study in Radiative Transfer has heen applied for some time to wa 
ter optics, by a large number researchers. The references cited are 
excellent reading for further information 

‘As mentioned, there is one additional radiosity scenario that can 
bbe important to ocean rendering under certain circumstances, but 
‘hich we will net consider. The stuaton stated gure 
Following the trail ofthe arrows, which track the direction light 
is travelling, we see that sometimes light coming to the surface 
(from above or below), ean reflect and/or transmit through the sur- 
face more than once. ‘The conditions which produce this behavior 
im significant amounts are: the wave heights must be fairly high, 
and the direction of viewing the waves, or the direction ofthe light 
Source must be nearly grazing the surface. ‘The higher the waves 
are, the less grazing the light source or camera need to be. This 
phenonmenon has been examined experimentally and in computer 
simulations. It is reasonably well understood, and we will ignore it 
from this point on 


3 Mathematics, Physics, and Experi- 
ments on the Motion for the Surface 


In this section we take a look at the mathematical problem we are 
tuying to solve. We simplify the mathematics considerably by ap- 
plying a series of approximations, How do we know these approx- 
mations are any good? There are decades of oceanographic re- 
search in which the ocean surface motion has been characterized 
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by measurements, simulations, mathematical analysis, and experi- 
mentation, The approximations we apply in this section are not per- 
fect, and there are many circumstances in the real world in which 
they break down, But they work extraordinarily well for most con- 
ditions at sea. To give you some idea of how well they work, we 
show some experimental work in sectionfi]that has been done which 
clearly shows these approximations at work in the real world. 


3.1 Bernoulli's Equation 


‘The starting point ofthe mathematical formulation of ocean surface 
motion is the incompressible Navier-Stokes equations. Chapter 2 
of (13) provides a thorough derivation of Bernoulli's equation from 
Navier-Stokes. We provide here the short version, and refer you to 
Kinsman’s or some other textbook for details. 

‘The incompressible Navier-Stokes equations for the velocity 
‘u(x, 1) ofa fluid a any postion x at any time tare 


Syuvu = —vp+F 0) 
Veu(x,t) = 0 (10) 


In these equations, p(x, t) is the pressure on the fluid, and F(x, t) 
is the force applied to the fluid. In our application, the force is 
conservative, so F = —VU for some potential energy function 
U(x,t), 

For ocean surface dynamics due to conservative forces like grav- 
ity, it turns out to be worthwhile to restrict the type of motion to 
class called potential flow. This is a situation in which the velocity 
has the form of a gradient: 


u(x,t) = V dx.) ay 


‘This restriction has the important effect of reducing the number de- 
grees of freedom of the flow from the three components of velocity 
to the single one of the potential velocity function @. In fact, us- 
ing this gradient form, the four Navier-Stokes equations transform 
into the two equations 


ae 
mR 


+h (voy = -P- (12) 


V' d(x,t) = 0 (3) 


Equation [Tis called Bernoulli's equation. As a fully nonlinear 
reduction of the Navier-Stokes equations, Bernoulli's equation is 
capable of simulating a variety of surface dynamics effects, inclad- 
ing wave breaking in shoaling shallow water (the bottom rises from 
deep water up toa beach). For more detail on numerical simulations 
of Bernoulli's equation in 3D, see (23). 


3.2 Linearization 


For our purposes, we want to reduce the complexity of Bernoulli's 
equation even further by applying two restrictions: linearize the 
equations of motion, and limit evaluation of the equations to just 
points on the surface itself, ignoring the volume below the sur- 
face. This may seem like an extreme restriction, but when com- 
bined with some phenomenolical knowledge of the ocean, this re- 
strictions work very wel. 

The first restriction is to linearize Bernoulli's equation, This is 
simply the task of removing the quadratic term 1/2(V6)*. Elimi- 
nating this term means that we are most likely restricted to surface 
‘waves that are not extremely violent in their motion, at least in prin- 
ciple, So Bernoulli's equation is reduced further to 


a9 


«sy 


All of the quantities ¢, p, and (are still evaluated at 3D points x 
fon the surface and in the water volume. 

‘The second restriction is to evaluation quantities only on the wa~ 
ter surface. To do this we have to first characterize what we mean by 
the surface. We will take the surface to be a dynamically changing 
height field, (x,t) that isa function of only the horizontal posi- 
tion xc. and time i. For convenience. we define the mean height of 
the wate surface as the zero value of the height. With this definition 
of wave height, the gravity-induced potential energy term U7 is 


U=gh as) 


and q is the gravity constant, usually 9.8 m/sec? in metric units. 

Restricting to just the water surface has several important conse- 
‘quences. One of the first consequences is for mass conservation. In 
the incompressible Navier-Stokes equation, mass is conserved via 
the mass flux equation 


V- u(x,t) =0 (16) 


‘When we chose to consider only potential low, this mass conserva- 
tion equation became 


V? o(x,t) =0 an 


If we label the horizontal portion of the position vector as x. , so 
that 2 = (2c1.,y), and y is the coordinate pointing down into the 
‘water volume, then the mass conservation equation restricted to the 
surface looks in more detail like 


{vt ie } oc, =0 Ct) 


Now, when you look at this equation and see that ¢ now depends 
only on the x, on the surface, you might be tempted to throw out 
the 0? /y? part of the equation, because there does not appear to 
bbe a dependence. That would produce useless results. Instead, what 
works better in this odd world of partial differential equations is to 
allow 6 to be an arbitrary function (at least with respect to this mass 
conservation equation) and to define the y-derivative operator to be 


4V-0E ) 


so that the operator V? is zero, We will use this approach for any 
‘quantity evaluated on the water surface whenever we need a ver 
cal derivative. Of course, this introduces an unusual operator that 
contains a square root function, 

‘Another consequence of restricting ourselves to just the surface 
is that the pressure remains essentially constant, and we can choose 
to have that constant be 0. With this and the rest of the restrictions, 
Bernoulli's equation has been linearized to 


a 
By 


dota.t) 
Ot 


= W(x.) 20) 


‘There is one final equation that must be rewritten for ths sita- 
ation. Recall thatthe velocity potential ¢ is used to compute the 
3D fluid velocity as a gradient, u = Vo. The vertical component 
of the velocity must now use equation(T9] In addition, the vertical 
velocity ofthe fluid is the sime at the speed ofthe surface height. 
Combining these we get 


AML) ST Kix e 
ee) = VT ole t) ay 


‘These last two equations, ZO]and [IT] are the final equations of 
motion that are needed to solve for the surface motion, They can 
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also be converted into a single equation. For example, if we take a 
derivative with respect to time of equation [ZT] and use[0]to sub- 
stitute for the time derivative of the velocity potential, we get the 
single equation for the evolution of the surface height. 


2) 


‘This still involves the unusual operator \/—=V7. However, tak- 
ing two more time derivatives converts it to a more normal two- 
dimensional Laplacian for the equation 


thief) _ 
oF 


4 hocs,t) 23) 


This form is frequently the starting point for building mathematical 
solutions to the surface wave equation, 


3.3 Dispersion 


1k turns out that the equations we have built for surface height 
‘whether in the form of equations/20]andP27] or equation[33] or equa- 
tio] reduce to one primary lesson about surface wave propga- 
tion, This lesson is embodied in a simple mathematical relationship 
called the Dispersion Relation, which is the focus of this section. 
Our goal here is to obtain that simple expression from the math- 
ematics above, understand some of its meaning, and demonstrate 
that, even though it appears unrealistically simplistic, the Disper- 
sion Relation is in fact present in natural ocean waves and can be 
measured experimentally 

Lets just use the version of the surface evolution equation in 
equation B3]for convenience. The other two versions could be used 
and arrive at the same answer using a slightly different set of ma- 
nipulations, Note that the equation of motion for the surface height 
is linear in the surface height. So as with any linear differential 
equation, the general solution of the equation is obtained by adding 
up any numberof specific solutions, Soles find a specfie solution, 
W turns out that all specifi solutions have the form 


A(x, f) = hg exp {ike-2, — iwt} 4) 


The 2D vector kcand the numbers w and ho are generic parameter 
at this point. If we use this form of a solution in equation [23} it 
tums into an algebraic equation ike this: 


hy {ut —@} =0 5) 


(and i is the magnitude of the vector kc). For this solution, there 
are only two possibilities: 


1. hy = 0. Then the surface height is flat, and the solution is not 
very interesting. 


£.V@F and ho can be anything, This isthe interesting 
solution. 


‘What we have found here is that the entire Navier-Stokes fluid d 
namics problem, reduced to an evolution equation for the water sur- 
face and approximated to something that can be solved relatively 
easily, amounts toa single equation imposing the constraint that the 
temporal frequency w of surface height movements is connected to 
the spatial extent ofa the propagating wave k = |k). This relation- 
ship, w = +VgF is the Dispersion Relation mentioned earlier. 
Note in particular that there is no constraint placed on the ampli- 
tude fig. But if the Navier-Stokes equation does not have anything 
to say about the amplitude, how do we give it a value? One way 
is by imposing initial conditions on the height and on its vertical 
speed. For ocean surface simulation in the next section, we will use 


Figure 3: A slice through a 3D PSD showing that the observed 
‘wave energy follows the deep water dispersion relation very well. 


4am altemate method, a statistical procedure to generate random real- 
izations of the amplitude, guided by measurements of the variance 
properties of wave height on the open ocean. 

So the question remaining is just how reasonable is the disper- 
sion relation for modelling realistic ocean surface waves? This is 
where lots of experimental research can come into play. Although 
there have been many decades of research on ocean wave proper- 
ties using devices placed in the water to directly measure the wave 
motion at a point, here we look at some relatively new research that 
involves measuring wave properties remotely with a camera in a 
plane. 

‘The AROSS (24} is a panchromatic camera mounted in a special 
hosing on the nose of a small airplane, Attached to the camera is 
navigation and GPS instrumentation which allow the camera po- 
sition, viewing direction, and orientation to be measured for each 
frame. After the plane flys a circular orbit around a spot over the 
‘ocean, this data can be used to remap images of the ocean into a 
‘common reference frame, so that the motion of the aireraft has been 
removed (except for lighting variations). This remapping allows 
the researchers to use many frames of ocean imagery, typically 1-2 
minutes worth, in some data processing to look for the dispersion 
relation, 

‘The data processing that AROSS imagery is subjected to gen- 
crates something called a 3D Power Spectral Density (PSD). This 
is obtained by taking the Fourier Transform of a time series of im- 
agery in time, as well as Fourier Transforms in the two spatial di 
rections of images. The output of these 3 Fourier Transforms is a 
quantity that is closely related to the amplitudes fi (Ke, w) for each 
spatial and temporal frequency. These are then absolute squared 
‘and smoothed or averaged in some way so that the output is a nu- 
merical approximation of a statistical average of |fg(K, w)|?. 

But how does a 3D PSD help us decide whether the dispersion 
relation appears in nature? If the imagery found only dispersion 
constrained surface waves, then the 3D PSD should have the value 
O for all values of k, w that do not satisfy the dispersion relation. So 
mostly we would expect the 3D PSD to only have significant values 
ina narrow set of k, wv values 

Figure[§]show a plot of the 3D PSD generated from AROSS im- 
ages (24). From the 3D volumetric PSD. this plot figure is a plane 
sliced through the volume. When sliced like this, the dispersion 
relation is a curve on the slice, shown as two dotted curves, ‘The 
data is plotted as contours of PSD intensity, color-coded by the key 
‘on the right. PSD levels following the dispersion relation curves 
are much higher than in other regions. This shows that the motion 
of the surface waves on all scales includes a very strong dispersion 
relation style of motion. There are other types of motion certainly. 
which the PSD figure shows as intensity levels away from the dis- 


4. PRACTICAL OCEAN WAVE ALGORITHMS 


Figure 4: Site at which video data was collected in 1986, near Zuma 
Beach, California, 


persion curves, But the dispersion motion is the strongest feature 
of this data, 

Relatively simple experiments can be done by anyone with ac- 
cess to a video camera and a hilltop overlook of an ocean, For 
example, figure[ljis a frame from a video segment showing wa- 
ter coming into the beach near Zuma Beach, California. The video 
camera was located on hill overlooking the beach, in 1986, In 1993, 
the region of video frames indicated in the figure was digitized, to 
produce a time series of frames containing just water surface. 

Figure[5]shows the actual 3D PSD from the image data. There 
are two clear branches along the dispersion relationship we have 
discussed, with no apparent modification by shallow water affects. 
There is also a third branch that is approximately a straight line 
lying between the first two. Examination of the video shows that 
this branch comes from a surfactant layer floating on the water in 
part of the video frame, and moving with a constant speed. Exclud- 
ing the surface layer, this data clearly demonstrates the validity of 
the dispersion relationship, and demonstrates the usefulness of the 
linearized model of surface waves. 


4 Practical Ocean Wave Algorithms 


In this section we focus on algorithms and practical steps to build- 
ing height fields for ocean waves. Although we will be occupied 
mostly by a method based on Fast Fourier Transforms (FFTs), we 
begin by introducing a simpler description called Gerstner Waves. 
This is a good starting point for several reasons: the mathematics is 
relatively light compared to FFTS, several important oceanographic 
concepts can be introduced, and they give us a chance to discuss 
wave animation. After this discussion of Gerstner waves, we go 
after the more complex FFT method, which produces wave height 
fields that are more realistic. These waves, called “linear waves” or 
“gravity waves" are a fairly realistic representation of typical waves 
oon the ocean when the weather is not too stormy. Linear waves are 
certainly not the whole story, and so we discuss also some meth- 
ods by which oceanographers expand the description to “nonlinear 
waves”, waves passing over a shallow bottom, and very tiny waves 
about one millimeter across called capillary waves. 

In the course of this discussion, we will see how quantities like 
windspeed, surface tension, and gravitational acceleration come 
into the practical implementation of the algorithms. 
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Figure 5: Slice from a 3D Power Spectral Density grayscale plot, 
from processed video data, 


4.1 Gerstner Waves 


Gerstner waves were frst found as an approximate solution to the 
fluid dynamic equations almost 200 years ago. There first appli- 
cation in computer graphics seems to be the work by Fournier and 
Reeves in 1986 (cited previously). The physical model is to de- 
scribe the surface in terms ofthe motion of individual points on the 
surface. Toa good approximation, points on the surface ofthe water 
{g0 through a circular motion as a wave passes by. IF point on the 
Lndisturbed surface is labelled xo = (0,20) and the undisturbed 


height is yo = 0, then as a single wave with amplitude A passes by, 
the point on the surface is displaced at time t to 
x = x0 (k/k)Asin(le- x0 — wt) 26) 
y= Acos(ke-x0 —wt) en 


In these expressions, the vector k, called the wavevector, is a 
horizontal vector that points in the direction of travel of the wave, 
‘and has magnitude f related to the length of the wave (A) by. 


|. (28) 


‘The frequency 1 is related to the wavevector, as discussed later, 

Figurefp]shows two example wave profiles, each with a different 
‘value of the dimensionless amplitude kA. For values kA < 1, the 
"wave is periodic and shows a steepening at the tops of the waves as 
kA approaches |. For kA > 1,a loop forms atthe tops of the wave, 
‘and the “insides of the wave surface are outside", not a particularly 
desirable or realistic effect. 

‘As presented so far, Gerstner waves are rather limited because 
they are a single sine wave horizontally and vertically. However, 
this can be generalized to a more complex profile by summing a 
set of sine waves. One picks a set of wavevectors ki, amplitudes 
Ai, frequencies ws, and phases 6., for i = 1,...,.V. to get the 
expressions 


= Xo — SO (e/h)Ay sin(le, x9 — wit +44) 29) 
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Figure 6; Profiles of two single-mode Gerstner waves, with differ- 
cnt relative amplitudes and wavelengths. 


| cos(ki “2x0 — wit + ds) 30) 


Figure[J]shows an example with three waves in the set. Interest- 
ing and complex shapes can be obtained in this way 


4.2. Animating Waves: The Dispersion Relation 


‘The animated behavior of Gerstner waves is determined by the set 
of frequencies w, chosen for each component. For water waves, 
there is a well-known relationship between these frequencies and 
the magnitude of their corresponding wavevectors, k,. In deep wa- 
tet, where the bottom may be ignored, that relationship is 


on 


The parameter g is the gravitational constant, nominally 
9.8m/sec?. This dispersion relationship holds for Gerstner waves, 
and also for the FFT-based waves introduced next. 

‘There are several conditions in which the dispersion relationship 
is modified. When the bottom is relatively shallow compared to 
the length of the waves. the bottom has a retarding affect on the 
waves. For a bottom at a depth D below the mean water level, the 
dispersion relation is 


k) = ghtanh(kD) (32) 


Notice that if the bottom is very deep, the behavior of the tanh 
function reduces this dispersion relation to the previous one. 

‘A second situation which modifies the dispersion relation is sur- 
face tension. Very small waves, with a wavelength of about | em or 
less, have an additional term: 


w*(k) = gh(1 + ML) (33) 


and the parameter L has units of length. Its magnitude is the scale 
for the surface tension to have effect. 

Using these dispersion relationships, it is very difficult to create 
‘a sequence of frames of water surface which for a continuous loop. 


eg 


Figure 7: Profile of a 3-mode Gerstner wave, 


In order to have the sequence repeat after a certain amount of time 
T for example, it is necessary that all frequencies be multiples of 
the basic frequence 


G4 


T 


However, when the wavevectors k are distributed on a regular lat- 
tice, itis impossible to arrange the dispersion-generated frequencies 
to also be on a uniform lattce with spacing ws, 

‘The solution to that is to not use the dispersion frequences, but 
instead a set that is close to them, For a given wavenumber k, we 


teethe Requeny 
ah) = [[ alee Gs) 


‘where [fa] means take the integer part ofthe value ofa, and w(K) is 
ny dispersion relationship of interest. The frequencies 2(h) are a 
‘quantization ofthe dispersion surface, and the animation ofthe wa- 
ter surface loops after atime T because the quantized frequencies 
are all integer multiples of wo. Figure[B]plots the original disper- 
sion curve, along with quantized dispersion curves fr two choices 
ofthe repeat time T. 


4.3 Statistical Wave Models and the Fourier Trans- 
form 


‘Oceanographic literature tends to downplay Gerstner waves asa re- 
alistic model of the ocean. Instead, statistical models are used, in 
combination with experimental observations. In the statistical mod- 
els, the wave height is considered a random variable of horizontal 
position and time, h(x, ) 

Statistical models are also based on the ability to decompose 
the wave height field as a sum of sine and cosine waves. The 
value of this decomposition is that the amplitudes of the waves 
have nice mathematical and statistical properties, making it sim- 
plerto build models. Computationally, the decomposition uses Fast 
Fourier Transforms (fits), which are a rapid method of evaluating 
the sums 
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Figure 8: A comparison of the continuous dispersion curve w 
VF and quantized dispersion curves, for repeat times of 20 sec- 
nds and 100 seconds. Note that fora longer repeat time, the quan- 
tized is a closer approximation to the original curve. 


‘The fft-based representation of a wave height field expresses the 
wave height fi(x, t) at the horizontal position x = (.r, ) as the sum 
of sinusoids with complex, time-dependent amplitudes: 


h(x. t) = 0 h(t) exp (ak x) (36) 


where ¢ is the time and Ic is a two-dimensional vector with com- 
ponents k = (ke,k-). ke = 2nn/Ly, k. = 2nm/L., and 
n and m are integers with bounds —N/2 <n < N/2 and 
“M/2 < m < M/2. ‘The fit process generates the height field 
at discrete points x — (nL, The value at other 
Points can also be obtained by switching toa discrete fourier trans- 
form, but under many circumstances this is unnecessary and is not 
applied here. ‘The height amplitude Fourier components, f(t). 
determine the structure ofthe surface. The remainder of ths sub- 
section is concemed with generating random sets of amplitudes in 
away that is consistent with oceanographic phenomenology. 

For computer graphics purposes, the slope yector of the wave- 
height field is also needed in order o find the surface normal. angles 
of incidence, and other aspects of optical modeling as weil. One 
‘way to compute the slope is though a finite difference between fit 
grid points, separated horizontally by some 2D vector Ax. While 
4 finite difference is efficient in terms of memory requirements, it 
can be a poor approximation to the slope of waves with small wave- 
length. An exact computation ofthe slope vector can be obtained 
by using more fis 


c(t) = Vib, 


ve Ak, t) exp (ike +x) (7) 


In terms of this fft representation, the finite difference approach 
‘would replace the term ik with terms proportional to 


exp (ike Ax) — 1 8) 


‘hich, for small wavelength waves, does not well approximate the 
gradient ofthe wave height. Whenever possible, slope computation 
Via the fit in equation[¥7}is the prefered method. 
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‘The fit representation produces waves on a patch with horizontal 
dimensions L., x L.., outside of which the surface is perfectly peri- 
‘dic. In practical applications, patch sizes vary from 10 meters to 2 
kilometers on a side, with the number of discrete sample points as 
hhigh as 2048 in each direction (i.e. grids that are 2048 x 2048, or 
‘over 4 million waves). The patch can be tiled seamlessly as desired 
over an area. The consequence of such atiled extension, however, is 
that an artificial periodicity in the wave field is present. As long as 
the patch size is large compared to the field of view, this periodicity 
is unnoticeable. Also, if the camera is near the surface so that the 
effective horizon is one or two patch lengths away, the periodicity 
‘will not be noticeable in the look-direction, but it may be apparent 
as repeated structures across the field of view. 

Oceanographic research has demonstrated that equation jis a 
reasonable representation of naturally occurring wind-waves in the 
‘open ocean. Statistical analysis of a number of wave-buoy, photo- 
‘graphic, and radar measurements of the ocean surface demonstrates 
that the wave height amplitudes /i(k, t) are nearly statistically sta- 
tionary, independent, gaussian fluctuations with a spatial spectrum 
denoted by 


x(k) = (|f*(,0)[*) 39) 


for data-estimated ensemble averages denoted by the brackets ( 
‘There are several analytical semi-empirical models forthe wave 

spectrum P;(k). A usefal model for wind-driven waves larger than 

capillary waves ina fully developed sea i the Phillips spectrum 


OO eat, ao 


Pals) = 


where L = V?/g is the largest possible waves arising from a con- 
tinuous wind of speed Vg is the gravitational constant, and 1 is 
the direction ofthe wind, isa numeric constant. The cosine factor 
[k-2i[? in the Phillips spectrum eliminates waves that move perpen- 
dicular to the wind direction. This model, while relatively simple, 
thas poor convergence properties at high values of the wavenumber 
||. A simple fix is to suppress waves smaller that a small length 
£& L, and modify the Phillips spectrum by the multiplicative fac- 
tor 


exp (-0) ay 
OF course, you are fee to “rll your own spectrum to try out 
various effect. 


4.4 Building a Random Ocean Wave Height Field 


Realizatons of water wave height lds ae created from the prin- 
ciples elaborated up to tis point) gaastan random nambers with 
Spatial spectra of a prescribe form. This is most eficiently accom 
lished directly inthe fourer domain. The fourier amplitudes of 
Grave height Sold cea be prodnced at 

otk) = + 6 


va i) VPi(k) (42) 


where &, and €, are ordinary independent draws from a gaussian 
random number generator, with mean O and standard deviation | 
Gaussian distributed random numbers tend to follow the experi- 
mental data on ocean waves, but of course other random number 
<istributions could be used. For example, log-normal distributions 
could he used to produce height fields that are vary “intermittent” 
i.e. the waves are very high or nearly fat, with relatively litle in 
between, 

Given a dispersion relation .o(), the Fourier amplitudes of the 
‘wave field realization at time ¢ are 


ho(k) exp {iw(k)t} 
+ fig(—k)exp{—in(k)t} 43) 


fut.) 
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‘This form preserves the complex conjugation property h” (kt) = 
h(t) by propagating waves “to the left” and “to the right”. In 
addition to being simple to implement, this expression is also effi- 
cient for computing f(x), since it relies on ffs, and because the 
‘wave field at any chosen time can be computed without computing 
the field at any other time. 

In practice, how big does the Fourier grid need to be? What 
range of scales is reasonable to choose? IF you want to generate 
wave heights faster, what do you do? Lets take a look at these 
questions. 


How big should the Fourier grid be? The values of N and M 
can be between 16 and 2048, in powers of two. For many 
situations, values in the range 128 to 512 are sufficient. For 
extremely detailed surfaces, 1024 and 2048 can be used. For 
example, the wave fields used in the motion pictures Water- 
world and Titanic were 2048x2048 in size, with the spacing 
between grid points at about 3.cm, Above a value of 2048, one 
should be careful because the limits of numerical accuracy for 
floating point calculations can become noticeable. 


What range of scales is reasonable to choose? The answer to this 
question comes down to choosing values for L, Lz, M, and 
NN, The smallest facet in either direction is dr = L,/M or 
Generally, dr and dz need never go below 
Below this scale, the amount of wave action is 
small compared to the rest of the waves. Also, the physics 
of wave behavior below 2 cm begins to take on a very differ- 
ent character, involving surface tension and “nonlinear” pro- 
cesses. From the form of the spectrum, waves with a wave- 
length larger than V2/g are suppressed. So make sure that dr 
and dz are smaller than V2/q by a substantial amount (10 ~ 
1000) or most of the interesting waves will be lost. The se- 
cret to realistic looking waves (c.g. figure[T2|(a) compared to 
figure[T2](c)) is to have MI and JV as large as reasonable. 


How do you generate wave height fields in the fastest time? The 
time consuming part of the computation is the fast fourier 
transform, Running on a 1+ GHz epu, 512 x 512 FFTs can 
be generated at nearly interactive rates, 


4.5 Examples: Height Fields and Renderings 


We now turn to some examples of waves created using the fft ap- 
proach discussed above, We will show waves in two formats: as 
greyscale images in which the grey level is proportional to wave 
height; and renderings of oceanscapes using several different ren- 
dering packages to illustrate what is possible. 

In the first set of examples, the grid size is set to M = N = 512, 
with Le 1000 meters. The wind speed is a gale force at 
V = 81 meters/second, moving in the x-direction. The small-wave 
cutoff of € = 1 meter was also used. Figure Dis a greyscale rep- 
resentation of the wave height: brighter means higher and darker 
means lower height. Although produced by the fft algorithms de- 
scribed here, figure[B]is not obviously a water height field. It may 
help to examine figure[TO] which is a greyscale depiction of the 
x-component of the slope. This looks more like water waves that 
figure[9] What is going on’? 

Figures| [demonstrate consequence of water surface 
optics, discussed im the next section: the visible qualities of the 
surface structure tend to be strongly influenced by the slope of the 
waves, We will discuss this in quantitative detail, but for now we 
willl summarize it by saying that the reflectivity of the water is a 
strong function of the slope of the waves, as well as the directions 
of the light(s) and camera. 

To illustrate a simple effect of customizing the spectrum model, 
figure[[T}is the greyscale display of a height field identical to figure 
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Figure 9: A surface wave height realization, displayed in greyscale, 


Figure 10: The x-component of the slope for the wave height real- 
ization in figure[] 
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Figure 11: Wave height realization with increased directional de: 
pendence. 


fp] with the exception that the directional factor {Ie -10[? in equation 
[o}has been changed to [K-1i|®, The surface is clearly more aligned 
‘With the direction of the wind. 

‘The next example of a height field uses a relatively simple shader 
in BMRT, the Renderman-compliant raytracer, The shader is shown 
in the next section. Figure I2]shows three renderings of water sur- 
faces, varying the size of the grid numbers M and N’ and making 
the fucet Sizes dx and dz proportional to 1/M and 1/N. So as we 
go from the top image to the bottom, the facet sizes become smaller, 
and we see the effect of increasing amount of detail in the render: 
ings. Clearly, more wave detail helps to build a realistic-looking 
surface. 

As a final example, figure[T3]is an image rendered in the com: 
mercial package RenderWorld by Arete Entertainment. This ren: 
dering includes the effect of an atmosphere, and water volume scat 
tered light. These are discussed in the next section. But clearly, 
wave height fields generated from random numbers using an fit pre 


scription can produce some nice images. 


4.6 Choppy Waves 


We tum briefly in this section to the subject of creating choppy 
looking waves. The waves produced by the fit methods presented 
up to this point have rounded peaks and troughs that give them the 
appearance of fair-weather conditions. Even in fairly good weather, 
and particularly in a good wind or storm, the waves are sharply 
peaked at their tops, and flattened at the bottoms. The extent of this 
chopping of the wave profile depends on the environmental condi 

tions, the wavelengths and heights of the waves. Waves that are 
sufficiently high (e.g. with a slope greater than about 1/6) eventu: 

ally break at the top. generating a new set of physical phenonema 
in foam, splash, bubbles, and spray. 

‘The starting point for this method is the fundamental fluid dy- 
namic equations of motion for the surface. These equations are ex. 
pressed in terms of two dynamical fields: the surface elevation and 
the velocity potential on the surface, and derive from the Navier: 
Stokes description of the fluid throughout the volume of the water 
and air, including both above and below the interface. Creamer 


ie 12: Rendering of waves with (top) a fairly low number of 
waves (facet size = 10 em), with little detail: (middle) a reasonably 
good number of waves (facet size = 5 em); (bottom) a high number 
‘of waves with the most detail (facet size = 2.5 cm), 


Figure 13: An image of a wave height field rendered in a commer. 
cial package with a model atmosphere and sophisticated shading, 
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Figure 14: A comparison of a wave height profile with and without 
the displacement. The dashed curve isthe wave height produced by 
the ft representation. The solid curve isthe height field displaced 
using equation] 


et alf{) set out to apply a mathematical approach called the "Lie 
Transform technique" to generate a sequence of “canonical trans- 
formations” of the elevation and velocity potential. The benefit of 
this complex mathematical procedure is o convert the elevation and 
velocity potential into new dynamical fields that have a simpler 
dynamics, The transformed case is in fact just the simple ocean 
height field we have been discussing, including evolution with the 
same dispersion relation we have been using in this paper. Start- 
ing from there, the inverse Lie Transform in principle converts our 
phenomenological solution into a dynamically more accurate one 
However. the Lie Transform is difficult to manipulate in 3 dimen- 
sions, while in two dimensions exact results have been obtained. 
Based on those exact results in two dimensions, an extrapolation 
for the form of the 3D solution has been proposed: a horizontal dis- 
placement of the waves, with the displacement locally varying with 
the waves. 

In the ft representation, the 2D displacement vector field is com- 
puted using the Fourier amplitudes of the height field, as 


peo-> . i(k, t) exp (ike-x) (44) 


Using this vector field, the horizontal position of a grid point of 
the surface is now x + AD(x,f), with height f(x) as before. The 
parameter A is not part of the original conjecture. but is a conve- 
nient method of scaling the importance of the displacement vector. 
This conjectured solution does not alter the wave heights directly, 
but instead warps the horizontal positions of the surface points 
a way that depends on the spatial structure of the height field. The 
particular form of this warping however, actually sharpens peaks in 
the height field and broadens valleys, which is the kind of nonlin- 
cear behavior that should make the fft representation more realistic. 
Figure[Tq] shows a profile of the wave height along one direction 
in a simulated surface. This clearly shows that the “displacement 
conjecture” can dramatically alter the surface. 

‘The displacment form of the this solution is similar to the algo- 
rithm for building Gerstner waves [I2} discussed in section] In 


Figure 15: A time sequence of profiles of a wave surface. From top 
to bottom, the time between profiles is 0.5 seconds. 


that case however, the displacement behavior, applied to sinusoid 
shapes, was the principle method of characterizing the water sur- 
face structure, and here it is a modifier to an already useful wave 
height representation. 

Frgure [5]illstes how these choppy waves behave a they 
evolve. The taps of waves form a sharp cusp, which rounds out 
and disappears shortly afterward, 

One “problem” with this method of generating choppy waves 
can be sen in gue [T] Nea the tops of some ofthe waves the 
surface actally passes Through itself and inverts, so thatthe outward 
normal tothe surface points inward. This is because the amplitudes 
of the wave components can be large enough to create large dis- 
placements that overlap. This is easily defeated simply by reducing 
the magnitude of the scaling factor \. For the purposes of computer 
‘graphics, this might actually be a useful effect to signal the pro- 
duction of spray, foam and/or breaking waves. We will not discuss 
hhere how to carry out such an extension, except to note that in order 
to use this region of overlap, a simple and quick test is needed for 
ddciding thatthe effect is taking place. Fortunately, there is such a 
simple test in the form of the Jacobian of the transformation from 
xox + AD(2x, t). The Jacobian is a measure of the uniqueness of 
the transformation. When the displacement is zero, the Jacobian is 
1. When there is displacement, the Jacobian has the form 


IK) = Ses Iyy = Sevdys 5 (45) 


with individual terms 


1g youl) 
Jee(x) = 14553 
Iyy(x) = 14 ,2Deb) 
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Figure 16: Wave height profile with and without the displacement, 
Aso plotted is the Jacobian map for choppy wave profile. 


Jay (X) 


and D = (Dz, D,)). The Jacobian signals the presence of the over- 
lapping wave bacause its value is less than zero in the overlap re- 
gion, For example, figure[[6)plots a profile of a basic wave without 
displacement, the wave with displacement, and the value of J for 
the choppy wave (labeled "Folding Map"). The “folds” or overlaps 
in the choppy surface are clearly visible, and align with the regions 
in which J < 0. With this information, it should be relatively easy 
to extract the overlapping region and use it for other purposes if 
desired. 

But there is more that can be learned from these folded waves 
from a closer examination of this folding criterion, The Jacobian 
derives from a2 x 2 matrix which measures the local uniqueness 
of the choppy wave map x — x + AD. This matrix can in general 
be written in terms of eigenvectors and eigenvalues as: 


Jou = Je + F,eref , (ab=2,y) (46) 


where .J- and J; are the two eigenvalues of the matrix, ordered 


so that J_ < J. The corresponding orthonormal eigenvectors are 
€~ and &* respectively. From this expression, the Jacobian is just 
Ja Ids. 


‘The criterion for folding that J < 0 means that J. < 0 and 
J > 0. So the minimum eigenvalue isthe actual signal of the on- 
Set of folding. Further, the eigenvector & points inthe horizontal 
direction in which the folding is taking place. So, the prescrip- 
tion now is to watch the minimum eigenvalue for when it becomes 
negative, and the alignment of the folded wave is parallel to the 
minimum eigenvector. 

‘We can illustrate this phenomenon with an example. Figures{T7] 
and[T¥]show two images of an ocean surface, one without choppy 
‘waves, and the other withthe choppy waves strongly applied. These 
two surfuces are identical except for the choppy wave algorithm. 
Figure T9]shows the wave profiles of both surfaces along a slice 
though the surfaces. Finally, the profile of the choppy wave is 
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plotted together with the value of the minimum eigenvalue in figure 
(20) showing the cen eommecton betwen folding and the negative 
value of J 

Incidentally, computing the eigenvalues and eigenvectors of this 
matrix is fast because they have analytic expressions as 


1 21/2 
wv) + 5 { (ex — Jun)? + 4ey} (47) 
for the eigenvalues and 

(as) 


with 
(49) 


for the eigenvectors, 


5 Interactive Waves from Disturbances 


‘The Fourier based approach to water surface evolution described in 
the section ]has several limitations that make in unworkable for 
really interactive applications, Fundamentally, the Fourier method 
computes the wave height everywhere in one FFT computation, 
You cannot choose to compute wave height in limited areas of a 
surface grid without completely altering the calculation, You either 
get the whole surface, or you don't compute it. This makes it very 
hard to customize the wave propagation problem from one location 
to another. So it you want to put some arbitrarily-shaped object 
in the water, move it around some user-constructed path, the FFT 
method makes it dificult to compute the wave response tothe shape 
and movement of the object. So if you want to have an odd looking 
craft moving around in the water, and maybe going up and down 
and changing shape, the FFT method is not the easiest thing to use. 
‘Also, if you want to have a shallow bottom, with a sloping beach or 
‘underwater sea mound or some other variability in the depth of the 
water, the FFT method again is a challenge use 

Fortunately in the last few years an alternative scheme has 
emerged which allows fast construction of a water surface in re- 
sponse to interactions with objects in the water and/or variable bot- 
tom depth, The heart of the method is an approach to computing 
the propagation which does not use the FFT method at al, The fun- 
damental mathematics of this interactive method is deseribed in the 
reference 22], and is refered to as ‘Wave, Here we very quickly run 
through the approach and show some examples in action 

‘Avits heart the iWave method returns tothe linearized equation 
for the wave height 22] We can turn the time derivatives imo finite 
time differences for sume step At, and get an explicit expression 
for the wave height 


h(x,t+At) = 2h(x,t)—h(x,t— At) 
= (AN? Y=VEA1) (50) 


‘This form can be used to explicitly advance the surface wave height 
from one frame to the next. The hard part of course i figuring out 
how to calculate the last term, with the square root of the Laplacian 
operator. 

‘The solution is to use convolution. The wave height is kept on 
a rectangular grid of points (he dimensions of which do not have 
to be powers of two), and we make use of the fact that any linear 
‘operation ona grid of data values can be converted into a convolM- 
tion of some sort. The details of the numerical implementation are 
contained in (23) 

But the real, amazing, property of iWave is that wave interaction 
with objects on the water surface is evaluated with a very simple 
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Figure 17: Simulated wave sur 
applied. Rendered in BMRT with 


ithout the choppy algorithm Figure 18: Same wave surface with strong chop applied. Rendered 
ic plastic shader. in BMRT will plastic shader. 
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Figure 19: Profiles of the two surfaces, showing the effect of the 
choppy mapping. 


Figure 20: Plot of the choppy surface profile and the minimum 
eigenvalue, The locations of folds of the surface are clearly the 
same as where the eigenvalue is negative, 
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5.1 Modifications for shallow water 


\Wave simulations based on the FFT method can simulate shallow 
water effects by using the dispersion relationship in equation [2] 
This only applies toa flar bottom. It would be nice to simulate wave 
propagation onto.a beach, or pasta shallow subsurface sea mount, oF 
‘vera submarine that i just below the surface. The iWave method 
is a great way of doing those things. Recall that the iWave method 
converts the mathematical operation 


oh 1) 


into a convolution, so that effectively g /=V? becomes a convolu- 
tion kernel. A shallow bottom with depth D changes this term to 
(compare with the dispersion relationship) 


1/77 tah(V=%D) bw 


This operation can also be converted into a convolution, with 
9V=V? tanh (V=V2D) becomes a convolution kernel. 

How is this applied to a variable-depth bottom? The convolution 
keel is a 13x13 matrix(®2). So we could build a collection of 
kemels over a range of depth values D. At each point on the grid, a 
custom kemel is constructed for the actual depth at that grid point 
by interpolating from the set of prebuilt values. Figure P3]is the 
‘wave height from a simulation using this technique. The bottom 
slopes from deep on the right hand sid, to a depth of O on the left 
edge. In addition, there isa subsurface sea mount on the right. 

‘There are three important behaviors inthis simulation that occur 
in real shallow water propagation: 


1 Waves in shallow regions have large amplitudes that in deep 


2. As waves approach a beach, they pile up together and have a 
higher spatial frequency. 


3. The subsurface sea mount causes a diffraction of waves, 


In addition to these capabilities, iWave can also compute other 
‘quantities, such as cuspy waves and surface velocity. 


5. INTERACTIVE WAVES FROM DISTURBANCES 


a) 


ib) 


Ke) 


Figure 21: A sequence of frames from an iWave simulation, show- 
ing waves reflecting off of objects in the water. The objects are the 
black regions. In the upper left there is also diffraction taking place 
as waves propagate the narrow channel and emerge in the corner 
region, 
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Figure 22: Frame from a simulation and rendering showing waves 
interacting with a ship. 


Figure 23: Frame from an iWave simulation with a variable depth 
shallow bottom. 
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6 Surface Wave Optics 


‘The optical behavior of the ocean surface is fairly well understood, 
at least for the kinds of quiescent wave structure that we consider 
in these notes. Fundamentally, the ocean surface is a near perfect 
specular reflector, with well-understand relectivity and transmis- 
sity functions. In this section these properties are summarized, and 
combined into a simple shader for Renderman. There are circum- 
stances when the surface does not appear to be a specular reflector. 
In particular, direct sunlight reflected from waves ata large distance 
from the camera appear to be spread out and made diffuse. This 
is due to the collection of waves that are smaller than the camera 
ccan resolve at large distances. The mechanism is somewhat similar 
to the underlying microscopic reflection mechnanisms in solid sur- 
faces that lead to the Torrance-Sparrow model of BRDFs. Although 
the study of glitter patterns in the ocean was pioneered by Cox and 
Munk many years ago, the first models of this BRDF behavior that 
Lam aware of were developed in the early 1980's. At the end ofthis 
section, we introduce the concepts and conditions, state the results, 
and ignore the in-between analysis and derivation. 

‘Throughout these notes, and particularly in this section, we 
nore one optical phenomenon completely: polarization. Polariza- 
tion effects can be strong at a boundary interface like a water sur- 
face. However, since most of computer graphics under considera- 
tion ignores polarization, we will continue in that long tradition. OF 
course, interested readers can find literature on polarization effects 
at the air-water interface. 


6.1 Specular Reflection and Transmission 


Rays of light incident from above or below at the air-water interface 
are split into two components: a transmitted ray continuing through 
the interface at a refracted angle, and a reflected ray. The intensity 
of each of these two rays is diminished by reflectivity and trans- 
missivity coefficients. Here we discussed the directions of the two 
outgoing rays. In the next subsection the coefficients are discussed, 


6.1.1 Reflection 


‘As is well known, in a perfect specular reflection the reflected ray 
and the incident ray have the same angle with respect tothe surface 
normal. This is true for all specular reflections (ignoring roughen- 
ing effects), regardless of the material. We build here a compact 
expression for the outgoing reflected ray. First, we need to build up 
some notation and geometric quantities 

The three-dimensional points on the ocean surface can be la- 
beled by the horizontal position and the waveheight h(x, ) as 


¥(x,t) =x + Fh(x,t), (53) 
where ¥ is the unit vector pointing straight up. At the point r, the 


normal to the surface is computed directly from the surface slope 
(x,t) = VACx,t) as 


(54) 


For a ray intersecting the surface at from direction fi, the diree- 
tion of the reflected ray can depend only on the incident direction 
and the surface normal. Also, as mentioned before, the angle be- 
tween the surface normal and the reflected ray must be the same 
as the angle between incident ray and the surface normal. You can 
verify for yourself thatthe reflected direction fy is 


(x,t) = fay — Bis(x, t) (fag (x,t) - fh) (35) 


Note that this expression is valid for incident ray directions on either 
side of the surface, 
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6.1.2 Transmission 


Unfortunately, the direction of the transmitted ray is not expressed 
as simply as for the reflected ray. In this case we have two guid- 
ing principles: the transmitted direction is dependent only on the 
surface normal and incident directions, and Snell's Law relating the 
sines of the angles of the incident and transmitted angles to the in- 
dices of refraction of the two materials. 

‘Suppose the incident ray is coming from one of the two media 
with index of refraction nj (for air, n = 1, for water, n = 4/3 
approximately), and the transmitted ray is in the medium with index 
of refraction 7. For the incident ray at angle @, to the normal, 


sin®, = /T— (iy Bis = | x As] (56) 


the transmitted ray will be at an angle @ with 
sind, = jit, x tis| on 
Snell's Law states that these two angles are related by 
ng sin 8; = ne sin (58) 
‘We now have all the pieces needed to derive the direction of trans- 
mission. The direction vector can only be a linear combination of 


1, and fis. It must satisfy Snell's Law, and it must be a unit vector 
(by definition). This is adequate to obtain the expression 


f(x, t) = Pea + Pst) ast, 1) (59) 
‘with the funtion T defined as 


T(x,t) 


+ (60) 


‘The plus sign is used in T’ when fi - fs <0, and the minus sign is 
used when iy »fis > 0 


6.2 Fresnel Reflect 


Accompanying the process of reflection and transmission through 
the interface is a pair of coefficients that describe their efficiency. 
‘The reflectivity Rand transmissivity Tare related by the constraint 
that no light is lost atthe interface. This leads to the relationship 


and Transmissivity 


R+T=1 (1) 


‘The derivation of the expressions for F and T is based on the elec- 
tromagnetic theory of dielectrics. We will not carry out the deriva- 
tions, but merely write down the solution 


Rléis,8,) = 2 {wt =) , wont} @) 


2 Sn@ +0) ” tan", 74) 


Figure[2a}s a plot ofthe reflectivity for rays of ight traveling down 
onto Water surface as a function ofthe angle of incidence to the 
surface, The plot extends from a grazing angle of 0 degrees to per- 
pendicular incidence at 90 degrees. As should be clear, variation of 
the reflectivity across an image is an important source of the “tex 
ture" of feel of water. Notice that reflectivity is «function of the 
angle of incidence relative to the wave normal, which in tur is di 

rectly related tothe slope of the surface. So we ean expect that a 
strong contributor to the texture of water surface is the pattern of 
slope, while variation ofthe wave height serves primarily asa wave 
hiding mechanism. ‘This is the quantitative explanation of why the 
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Figure 24: Reflectivity for light coming from the air down to the 
water surface, as a function of the angle of incidence of the light. 


surface slope more closely resembles rendered water than the wave 
height does, as we saw in the previous section when discussing fig- 
wre] 

‘When the incident ray comes from below the water surface, there 
are important differences in the reflectivity and transmissivity. Fig- 
ureB5]shows the reflectivity asa function of incidence angle again, 
but this time for incident light from below. In this ease, the re- 
flectivity reaches unity at a fairly large angle, near 41 degrees. At 
incidence angles below that, the reflectivity is one and so there is no 
transmission of light through the interface. This phenomenon is o- 
tal internal reflection, and can be seen just by swimming around in 
4 pool. The angle at which total internal reflection begins is called 
Brewster's angle, and is given by. from Snell's Law, 


-™ 0.75 (63) 


sino? 


or a? 


48.6log. In our plots this angle is 90 — 02 


ALLdeg, 


6.3 Building a Shader for Renderman 


From the discussion so far, one of the most important features 3 ren= 
dering must emulate is the textures of the surface due to the strong 
slope-dependence of reflectivity and transmissivity. In this section 
‘we construct a simple Renderman-compliant shader using just these 
features. Readers who have experience with shaders will know how 
to extend this one immediately. 

‘The shader exploits that fact that the Renderman interface al- 
ready provides a built-in Fresnel quantity calculator, which pro- 
vides F, T. tip, and thy using the surface normal, incident direction 
vector, and index of refraction. The shader for the air-to-water case 
is as follows: 


surface watercolorshader ( 

upwelling = color(0, 0.2, 0.3) 
sky color(0.69, 0-84, 1) 
air color (0.2, 0.1,0.1)5 


string envmap 
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Figure 25: Reflectivity for light coming from below the water sur- 
face, as a function of the angle of incidence of the light. 


ctivitys 
normalize(I) 
vector all = normalize (Ng); 

Float costhetai = abs(nl . aN); 
Float thetai = acos(costhatai) ; 
float sinthetat = sin(thetai) /nSnel. 
Float thetat = asin(sinthetat); 

if (thetai == 0.0) 


else 
C 
float fs (thetat - thetaiy 
(chetat + thetai); 
Float ts = tan(thetat - thetai) 
/ tan(thetat + thetai); 
reflectivity = 0.5 " ( fs¥fs + tarts}; 


) 
vector dPE = P-E; 

ength(dPe) * Kdiffuse; 
dist = exp(-dist); 


float dist = 


Ci = dist * ( reflectivity * sky 
+ (I-reflectivity) * upwelling ) 
+ (Gedist)* airy 


‘There are two contributions to the color: light coming downward 
‘onto the surface with the default color of the sky, and light coming 
upward from the depths with a default color. This second term will 
be discussed in the next section. Itis important for incidence angles 
that are high in the sky, because the reflectivity is low and transmis- 
sivity is high, 


7 WATER VOLUME EFFECTS 


‘This shader was used to render the image in figure 
BMT raytrace renderer. For reference, the exact sat 
been rendered in[P6] with a generic plastic shader. Note that the 
realistic water shader tends to highlight the tops of the the waves, 
where the angle of incidence is nearly 90 degrees grazing and the 
reflectivity is high, while the sides of the waves are dark, where 
angle of incidence is nearly 0 that the reflectivity is low: 


7 Water Volume Effects 


‘The previous section was devoted to a discussion of the optical be- 
havior of the surface of the acean. In this section we focus on the 
optical behavior of the water volume below the surface. We begin 
with a discussion of the major optical effects the water volume has 
on light, followed by an introduction to color models researchers 
have built to try to connect the ocean color on any given day to un- 
derlying biological and physical processes, These models are built 
upon many years of in-situ measures off of ships and peers. Fi- 
nally, we discuss two important effects, caustics and sunbeams, that 
sometimes are hard to grasp, and which produce beautiful images 
when properly simulated, 


7.1 Scattering, Transmission, and Reflection by 
the Water Volume 


In the open ocean, light is both scattering and absorhed by the vol- 
ume of the water, The sources for these events are of three types: 
water molecules, living and dead organic matter, and non-organic 
matter. In most oceans around the world, away from the shore lines, 
absorption is a fairly even mixture of water molecules and organic 
‘matter. Scattering is dominated by organic matter however. 

To simulate the processes of volumetric absorption and scatter- 
ing, there are five quantities that are of interest: absorption coef 
ficient, scattering coefficient, extinction coefficient, diffuse extinc- 
tion coefficient, and bulk reflectivity. All of these coefficients have 
units of inverse length, and represent the exponential rate of atten- 
uation of light with distance through the medium. The absorption 
coefficient « is the rate of absorption of light with distance, the 
scattering coefficient b is the rate of scattering with length, the ex- 
tinction coefficient cis the sum of the two previous ones c= a+b, 
and the diffuse extinction coefficient K° describes the rate of loss 
of intensity of light with distance after taking into account both ab- 
sorption and scattering processes. The connection between K and 
the other parameters is not completely understood, in part because 
there are a variety of ways to define Kin terms of operational mea- 
surements. Different ways change the details of the dependence. 
However, there is a condition called the asymptotic limit at very 
deep depths in the water, at which all operational definitions of KC 
converge to a single value. This asymptotic value of KC has been 
modeled in a variety of ways. There is a mathematically precise 
result that the ratio K’/e depends only on b/c, the single scatter 
albedo, and some details of the angular distribution of individual 
scattering distributions. Figure PJ8]is an example of a model of 
Ke for reasonable water condiffons. Models have been gener- 
ated for the color dependence of AC, most notably by Jerlov. In 
1990, Austin and Petzold performed a revised analysis of spectral 
models, including new data, to produce refined models of KC as a 
function of color. For typical visible light conditions in the ocean, 
K ranges in value from 0.03/meter to 0.L/meter. It is generally true 
thata < KC <c. 

‘One way to interpret these quantities for a simulation of water 
volume effects is as follows 


1. A ray of sunlight enters the water with intensity (after los- 
ing some intensity to Fresnel transmission). Along a path un- 
derwater of a length s, the intensity at the end of the path is 
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Figure 28: Dependence of the Diffuse Extinction Coefficient on the 
Single Scatter Albedo, normalized to the extinction, 


Texp(—es), ie. the ray of direct sunlight is attenuated as fast 
a possible 


2. Along the path through the water, a fraction of the ray is scat- 
tered into a distribution of directions. The strength of the scat- 
tering per unit length of the ray is 6, so the intensity is propor- 
tional to bf exp(—es). 


3. The light that is scattered out of the ray goes through poten- 
tially many more scattering events. It would be nearly im- 
possible to track all of them. However, the sum whole out- 
come of this process is to attenuate the ray along the path from 
the original path to the camera as bf exp(—cs) exp(—Kse), 
where sq is the distance from the scatter point in the ocean to 
the camera, 


A fifth quantity of interest for simulation is the bulk reflectivity 
of the water volume, This is a quantity that is intended to allow 
us to ignore the details of what is going on, treat the volume as a 
Lambertian reflector, and compute a value for bulk reflectivity. That 
number is sensitive to many factors, including wave surface condi- 
tions, sun angle, water optical properties, and details of the angular 
spread. Nevertheless, values of reflectivity around 0.04 seem to 
‘agree well with data, 


7.2 The Underwater POV: Refracted Skylight, 
Caustics, and Sunbeams 


Now that we have underwater optical properties at hand, we can 
look at two important phenomena in the ocean: caustics and sun- 
beams. 


7.2.1 Caustics 


Caustics, in this context, are a light pattern that is formed on sur- 
faces underwater, Because the water surface is not flat, groups of 
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Figure 29: Rendering of a caustic pattern at a shallow depth (5 
meters) below the surface. 


light rays incident on the surface are either focussed or defocussed 
‘Asa result, a point ona fictitious plane some depth below the ocean 
surface receives direct sunlight from several different positions on 
the surface. The intensity of light varies due to the depth, orig- 
inal contrast, and other factors. For now, lets write the intensity 
of the pattern as J = Ref Io, with Jy as the light intensity just 
above the water surface, The quantity Re isthe scaling factor that 
varies with position on the fictitious plane due to focussing and de- 
focussing of waves, and is called a caustic pattern, Figure3shows 
an example of the caustic pattern Ref. Notice thatthe caustic pat 
tem exhibits filaments and ring-like structure. Ata very deep depth 
the caustic patter is even more striking, as shown in figure 

‘One of the important properties of underwater light that produce 
caustic patterns is conservation of flux. This is actually a simple 
idea: suppose a small area on the ocean Surface has sunlight passing 
through it into the water, with intensity J at the surface. As we 
map that area to greater depths, the amount of area within it grows, 
or shrinks, but most likely grows depending on whether the area 
is focussed or defocussed. The intensity at any depth within the 
water is proportional to inverse of the area of the projected region 
Another way of saying ths is that if a bundle of light rays diverges, 
their intensities are reduced to keep the product of intensity time 
area fixed. 

Simulated caustic patterns can actually be compared (roughly) 
with real-world data, Ina series of papers published throughout the 
1970's, 1980's, and into the 1990's, Dera and others collected high 
speed time series of light intensity(2T]. As part ofthis data collec- 
tion and analysis project, the data was used to generate a probability 
distribution function (PDE) for the light intensity. Figure[ST]shows 
two PDFs taken from one of Dera's papers. The two PDF's were 
collected for different surface roughness conditions: rougher wa- 
ter tended to suppress more of the high magnitude fluctuations in 
intensity 

Figure[iJ}shows the pdf at two depths from a simulation of the 
because of many factors, but most importantly because we have 
not simulated the environmental conditions and instrumentation in 
Dera's experiments. Nevertheless, the similarity of figure[32}with 
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Figure 30; Rendering of a caustic pattern at great depth (100 me. 
ters) below the surface. 


Figure 31: PDF's as measured by Dera in reference (JT. 


Dera’s data is an encouraging point of information for the realism 
of the simulation. 


7.2.2 Godrays 


Underwater sunbeams, also called godrays, have avery similar or 
gin to caustics. Direct sunlight passes into the water volume, fo 
cussed and defocussed at different points across the surface. As the 
rays of light pass down through the volume. some of the light is 
scattered in other directions, and a fraction arives atthe camera 
‘The accumulated pattem of scatered light apparent to the camera 
are the godrays. So, while caustics are the pattem of direct sun 
light that penetrates down tothe floor of water volume, sunbeams 
are scattered light coming from those shafts of direct sunlight in 
the water volume, Figure[S3]demonstates sun bya 
camera looking up att 
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8 Appendix: Sample code for interactive 
water surfaces 


‘The code listed inthis appendix isa working implementation of the 
iWave algorithm. As listed below, iwvave paint looks like a crude 
paint routine. The user can paint in two modes. When iwave_paint 
Sars up. it begins running an iWave simulation on a small grid 
(200x200 withthe settings listed). This grid is small enough that 
the {Wave simulation runs interactively on cpus around 1 GHz. and 
faster. The running ofthe simulation is not apparent since there are 
no disturbing waves at start up. The interface looks like figure(34] 
In the default mode at start up, the painting generates obstructions 
that show up in black and block water waves, Figure [35] shows 
an example obstruction that has been painted, Hitting the 's’ key 
changes the painting mode to painting a source disturbance on the 
‘water surface, As you paint the disturbance, iwave paint propagates 
the disturbance, reflecting off of any obstructions you may have 
painted, Figure [36]shows a frame after source has been painted 
inside the obstructed and allows a brief period of time to propagate 
inside the obstruction and exit, creating a diffraction pattern at the 
mouth of the obstruction, 
‘A few useful keyboard options 


© Puts the paint mode into obstruction painting. This may be se- 
lected at any time. 


's Puts the paint mode into source painting. This may be selected 
at any time, 
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lave Demo 


Figure 35: The fwave paint window with obstruction painted. 


Wave Demo 


Figure 36: The ivave paint window after painting source inside the 


obstruction and letting it propagate out. 


Figure 34: The iwave.paint window at startup. 
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LA] = O.56(bmigh I /acal Sng tartar 41.0 [ebateuetbnt ily 


te Andee = a nah beget 


Mart betacttiniges 
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